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Abstract 

Nonlinear acoustics of wind instruments conducts to study unidimensional 
fluid flows. From physically relevant approximations that are modelized with the 
thin layer Navier Stokes equations, we propose a coupled model where perfect 
fluid flow is described by the Euler equations of gas dynamics and viscous and 
thermal boundary layer is modelized by a linear equation. We describe numeri- 
cal discretization, validate the associated software by comparison with analytical 
solutions and consider musical application of strongly nonlinear waves in the trom- 
bone. 
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Resume 

L'acoustique non lineaire des instruments a vents conduit a etudier les ecou- 
lements filaires monodimensionnels. A partir d'approximations physiquement 
realistes qui sont prises en compte par les equations de Navier Stokes de couche 
mince, nous proposons un modele couple ou le fluide parfait est deer it par les 
equations d'Euler de la dynamique des gaz et le fluide visqueux et conducteur de 
chaleur par une equation lineaire de couche limite. Nous detaillons la discretisation 
numerique retenue et validons le logiciel developpe grace a des solutions analy- 
tiques avant d'aborder I'application musicale a la propagation d'ondes fortement 
non lineaires dans le trombone. 



Key words : fluid mechanics, nonlinear acoustics, Euler equations, boundary 
layer, finite differences. 
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1) Introduction. 

• In this paper, we study simple models of nonlinear acoustic flows in cylindric or 
axisymmetric ducts. Our objective is to take into account several physical effects 
such compressibility of the air, viscous dissipation and thermal conduction, expe- 
cially in the vicinity of the wall. We consider the flow of a newtonian compressible 
fluid in a two-dimensional pipe. In a first approximation, the variation of physi- 
cal fields in the transverse direction can be neglected and an appropriate physical 
model for such a flow is given by unidimensional equations of gas dynamics (see, 
e.g. Landau-Lifschitz [LL53]). This model is appropriate for the description of 
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nonlinear waves in shock dynamics (see e.g. Courant-Friedrichs [CF48]) and also 
for weaker waves in nonlinear acoustics (Whitham [Wli74]). Nevertheless, such a 
model neglects all phenomena that can appear in the boundary layer. 

• The boundary layer is the region located near the wall (at a distance of the 
order of the boundary layer thickness 6) where viscous dissipation and thermal 
conduction have to be considered. The important role of the boundary layer in 
duct acoustics has been studied by Chester [Ch64]. Usually, linearized boundary 
layer equations are considered and the associated mathematical model is the heat 
equation whose solution can be explicited using a convolution kernel in time. 

• In the domain of acoustics, nonlinear and dissipative effects are usually taken 
into account via generalized Burgers equations as suggested by Blackstock [B185] ; 
this scalar model contains a source term which is, in the case of ducts, a convolution 
kernel giving an explicit solution of the linear model of acoustic boundary layer. 
We refer to Makarov-Ochmann [M097] for a review of the fundamental results. 

• We here focus on the fact that the modelling of an acoustic flow in a pipe can 
be conducted as a coupling between a perfect fluid and a boundary layer. We refer 
to Le Balleur [LB80] , Zeytounian [Ze92] and Aupoix-Brazier-Cousteix [ABC92] for 
classical approaches developed in the context of aerodynamics applications. In this 
paper, we have been inspired by these coupling techniques for pipe flow problem 
in nonlinear acoustics. 

2) Thin Layer Navier Stokes equations. 

• We consider the geometry of a two-dimensional pipe of characteristic longitu- 
dinal length equal to L. The thickness of the duct is 2/i and our basic hypothesis 

is that the ratio — is small : 
(2.1) J « 1. 

We consider also some longitudinal length A and some length / in the transverse 
direction {0 < I < h) that are used for the adimensionalization of the equa- 
tions (see Figure 1). We distinguish between the two types of flows associated to 
aerodynamics and acoustics applications. 
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Figure 1 Channel with characteristic lenghts. 

• In aerodynamics, we suppose simply : 

(2.2) A = L. 

Moreover, the distance I is a distance characteristic of the maximum of the bound- 
ary layer thickness (see e.g. Schlichting [Sc55]) : 

(2.3) l^,^ . 

In previous expression, p, fi and U are respectively the density, the viscosity 
and the amount velocity of the flow and we introduce also the so-called Reynolds 
number in aerodynamics : 

(2.4) 



We note that for extremely thin pipes, the boundary layer occupies all the duct, 
that is I c:^ h. In all cases, we suppose that 

2.5) e = - « , « 1 

A ^/W^ 



• In acoustics, if condition 

(2.6) ^i 

is satisfied, the waves propagate only along the axial direction (Pierce [Pi81], see 
also Bruneau [Br98]) and it is natural to consider the length wave A as a reference 
length for the axial direction. We set : 

(2.7) A = A. 

On the other hand, distance I is the natural length constructed from viscosity 
coefficient fi , density of the air pq and sound celerity cq at usual thermodynamic 
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conditions for pressure and temperature : po = 1 atmosphere and Tq = 300 
Kelvins. We set 

(2.8) / = 

po Co 

and we introduce also the acoustic Reynolds number R^^"^^ defined simirarily to 
expression (2.4) : 

(2.9) = 

The ratio — between right hand sides of expressions (2.8) and (2.7) satisfies the 
following hypothesis : 

/ 1 

(2.10 e = - « « 1. 

V / ^ ^acou 

In practice, I 10~^ m and if the frequency of the acoustic wave (with wave 
length A) is less than 1 Ghz, condition (2.10) is satisfied. More precisely, we set 
with Bruneau, Herzog, Kergomard and Polak [BHKP89], 

(2.11) Ivh = f^M+Mx;)— + (7-1) ^ 

where fiy is the volumic viscosity, 7 = 7/5 is the ratio of specific heats, k the 
thermic conductivity and Cp the calorific capacity at constant pressure. In the air 
the volumic viscosity fiy is negligeable compared to viscosity fi and (7 — 1) /c / Cp 
is of the order of viscosity /i i.e. the Prandtl number (see e.g. Schlichting [Sc55]) 
is of the order of 1 . Therefore, 

(2.12) I « U. 

that enforces hypothesis (2.8). 

• The fiow is supposed to satisfy the Navier Stokes equations of conservation 
of mass, impulse and energy. Recall that the unknowns are density velocity 
{u,v), pressure p and internal specific energy e. The thermodynamic variables 
are supposed to satisfy the state equation for the air that takes the classical form 
of a perfect gas equation 

(2.13) p = i7-l)pe. 

In the following, we neglect the volumic viscosity fiy and assume that the Stokes 
hypothesis concerning the two viscosities is valid. In consequence, the classical 
analytic expression of the Navier-Stokes equations {e.g. Landau and Lifschitz 
[LL53]) takes the form 
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(2.15) 
(2.16) 



d /4du 
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1 dv\ d"^ u 



dx\2> dx 3 dyJ dy 
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dx'^ ' dy^Sdx ' 3 dy 



+ 



(2.17) < 



_d_ 

dy 



+ pv ] = 



d_ 

dy 



= 1^ 

/ dv du^ 
\dx dy ^ 



d_ 

dx 



Adu 2 dv 



u - 



+ v\ - 



3 dx 3 dy 
Adv 2 du\ 



3 dy 3 dx 




• We detail the way we adimensionalize the Navier Stokes equations. First 
we have two length scales A and / for longitudinal and transverse directions 
respectively ; we denote by x and y these two space variables without dimension : 



(2.18) 
(2.19) 



_ X 

^ = A 
V 

y = r 



Second, we introduce some longitudinal reference velocity U . This velocity defines 
a time reference r and an adimensionalized time t according to 

A 

(2.20) 



(2.21) 



U 
t 



t = -. 

T 

If A = A is the length wave and U = cq is a typical choice for the adimension- 
nalization of velocity in acoustics, then r is the period of the wave, i.e. the time 
for the acoustic perturbation to travel one length wave. We introduce a second 
reference velocity V associated to this time r and the transverse distance / : 



(2.22) 



V = -. 

T 

If a particle travels distance A with axial velocity U during time r , it travels 
distance / with transverse velocity V during the same time interval. Therefore, 
we define dimensionless velocities u and v according to 



(2.23) 
(2.24) 



_ u 

" = u 



V 



V 

V 



1 V 
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with 

(2.25) e = |- . 

• For the adimensionahzation of convective terms, the reference for density is 
the density po of the air at the usual conditions and reference for pressure is 
associated with the dynamic pressure po U'^ . We set 

(2.26) p = — 

po 

^ = i^- 

The reference for internal energy is chosen in order to maintain the validity of the 
state equation (2.13) after adimensionahzation. We set 

(2.28) e = ^ 

and we deduce from (2.13), (2.26) (2.27) and (2.28) the state equation between 
these new variables : 

(2.29) p = (7-l)pe. 

• The Reynolds number Rq appears from the dissipation terms in the momen- 
tum equations (2.15) and (2.16) 

(2.30) fle = ^ . 

the Prandtl number P,. is defined from the heat fluxes in the energy equation 
(2.17) 



(2.31) Pr = 



k 

and a reference scale for temperature is defined by — : 

Cp 

(2.32) T = ^ . 

The Joule-Gay Lussac law for polytropic gas can be rewritten in terms of dimen- 
sionless energy e and temperature T according to 

(2.33) e = 7T. 

• Then the adimensionalized Navier Stokes equations take the following form 

(2.34) ^ + ^{pu) + {pv) = 
at ox ^ ' oy ^ ' 
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(2.35) 



(2.36) 



(2.37) < 



d d d 
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• The boundary conditions associated with these equations are of Dirichlet type 
on the boundary of the pipe : 



(2.38) 


u {x, 


y = 


-h) 


= u{x,y = 


h) = 


(2.39) 


v{x, 


y = 


-h) 


= v{x,y = 


h) = 


(2.40) 


T{x 


,2/ = 


-h) 


= Tix,y = 


--h) = To 



where To is the nondimensionless temperature given at the boundary of the pipe. 

• We observe first that the velocities u and v have the same order of magnitude 
and in consequence, due to the fact that 

(2.41) « 1 

we can neglect in the left hand side of equation (2.37) the v term compared to 
the u term. 



• We make the hypothesis that a typical distance for longitudinal variation of 

du U 

physical fields is of the order A . In particular — — ~ — and in consequence 

dx A 
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(2.42) g.l. 

We observe that no bigger gradients than — are taken in consideration into hy- 
pothesis (2.42) which means that the flow is regular and that no turbulence occurs. 
In an analogous way, a typical distance for transversal variation of all the fields is 

dv V U 

of the order of / and in particular 7— c^l — = — and we have again 

oy I A 

(2.43) 

dy 

d^ w 

More generally all the differential expressions of the type — 377 with w equal to 

d z'^ 

one of the physical fields p,u,v,T,e , variable z equal to t,x, or y and k — 1,2, 
is finally of the order of 1 (see e.g. Schlichting [Sc55] or Cousteix [C088]) : 

d'^ w 
'd¥ 



(2-44) ^ - 1 



• When we sum linear combinations of such expressions with coefficients of the 
type 1 or ^ (as in the right hand side of relation (2.35)), the leading term is the 
one that has the dominant factor ^ as a coefficient. We neglect in the following 
all the other terms. After these approximations, we have derived the so-called 
Thin Layer Navier Stokes equations (see e.g. Baldwin-Lomax [BL78], Kutler- 
Chakravarthy-Lombard [KCL78] or Rubin and Tannehill [RT92]). We re- write 
them without any adimensionalization : 

d d d d"^ u 

(2.46) —{pu) + ^(pu^ +p) + ^(puv) = fi 



dt^^ ' dx^^" ' ' dy^^^'' ^ dy"^ 

2.47 ^ = fi— -— + - — 

oy oy \ 3 ox 3 oy J 



(2.48) 



d ( / 1 2\ \ a / au\ , a^T 
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3) Perfect fluid for main flow. 

• We suppose now that the flow in the pipe satisfies the thin layer Navier-Stokes 
equations (2.45)-(2.48) and for fixed time t and abscissa x , we integrate equation 
(2.45) between y = —h and y = -\-h . We obtain in this way 



(3.1) 



i 



Due to boundary condition (2.39), the third term in (3.1) is null. We introduce 
now the mean values of density, momentum and energy in each x section according 
to : 

1 

Pit:X) = — J ^ p{t,x,y)dy 

1 

p{t,x)u{t,x) = ^{pu){t,x,y)dy 

1__ 1 1 

p{t, x) e{t, x) + -p{t, x) ^) = ^ / ( + -pw^ ) (t , x, y) dy . 

In terms of these new variables, the conservation of mass stands as : 



(3.2) 
(3.3) 
(3.4) 



(3.5) 



dp d 



at 







• In a similar way we integrate the impulse and the energy equations in the 
thickness of the pipe. We get 



(3.6) { 



d_ 
dt 



(3.7) I 



^(pe+ \pu') + j_^[pu{e+\)+vu){t^x,y)dy^ 



k f dT ^ 



dT 
dy 



{t,x,-h) 



due to boundary conditions (2.38) and (2.39). 



• We make first the hypothesis that the fields are quasi-constant in each section 
of the pipe and second that they have a rapid variation in a boundary layer region 
of thickness 6 , with the condition 
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(3.8) = « 1. 

The first hypothesis is absolutly non trivial. 

• In aerodynamics, it conducts (see e.g. Whitham [Wh74] and Msallam [Ms98]) 
to the shallow water equations when the following physical hypothesis is satisfied : 

(3-9) — — ^ — — = e— « 1 

due to hypothesis (2.1) and choice of variable e done in (2.5). 

• In acoustics, the Reynolds number is modified according to relation (2.9), i.e. 

(3.10) i?e-- = 

2 TT 

and the hypothesis of quasi-constancy of all the fields in the main flow is satisfied 
under the hypothesis (Kergomard [KeSl], Menguy and Gilbert [MG97]) 

1 1 

(3.11 — « 1. 

• This hypothesis can be justied as follow. Observe first that for a simple linear 
wave, the variation of pressure Pa due to acoustics satisfies the relation 

(3.12) Pa = PoCqu . 

dp 

Second the transverse gradient of pressure — satisfies at the first order a lin- 

oy 

earized version of equation (2.36) : 
(3.13) 



dp 




dv 


dy 







and the right hand side of this equation (3.13) can be evaluated as follow : 

dp 5 1 /(5\21 

(3.14) - « po-u- = pouco[j) J 

because r = — where the thickness of the boundary layer 8 is (see e.g. Lighthill 
Co ^ 

[Li78]) of the order of —= • 

(3.15) '5 ^ ^ 



acou 
e 



acou 
e 



We insert relations (3.8), (3.12) and (3.15) inside (3.14) and obtain 



(3.16) 
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dp Pa I I 



dy h Sh ' 

The transverse variations of pressure are of the order of the axial variations of 

pressure multiphed by the factor ^g^^ou g " Then relation (3.11) express that the 

transverse variation of pressure can be neglected compared with the axial ones. 

• We observe also that 

1 1 _ /Sy h _ 6h _ 6 h 1 1 

due to hypotheses (2.6), (3.8) and (3.15). Then hypothesis (3.11) is established 
even if the boundary layer thickness 6 is greater than the order of magnitude of 
the characteristic length I of visco-thermic effects. Physically, it corresponds to 
neglect volume losses compared to wall losses. 

• Under the hypothesis that all the fields are constant in the section 

-{h-S) < y < {h-S), 

we first observe that pressure p associated via the state equation (2.13) to mean 
density p and mean internal energy e can be well approached by the mean value 
of pressure : 

1 

(3.18) (7-l)pe ^ — / p{t,x,y)dy 

as mentioned previously. In an analogous way, we have 

1 

(3.19) pu^ ~ —J^^{pu^){t,x,y)dy 

(3.20) pu{e+-v?)+pu - ^ (pw(e+ y ) +pixj(t,a:,?/)d2/. 

All the hypotheses (3. 18)- (3.20) suppose finally that mean values of a nonlinear 
function is quasi-equal to the same nonlinear function of the mean values. This 
hypothesis is correct when the nonlinear function is well approximated by a con- 
stant. 

• We can now insert relations (3.18) to (3.20) inside equations (3.6) and (3.7). 
We obtain the final model for unidimensional perfect flow : 

(3.21) p(t^x) = (7-1)?? 
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(3.23) 



(3.24) 



11 f du , 



+ ^(^^'+^; 2h \dy 

-(pe+-pu ) + —[^pu{e+-u)+pu 

k f dT 



-(t,a;,-ft) 



dT 



4) Acoustic boundary layer. 

• We suppose as previously that the flow in the pipe satisfies the Thin Layer 
Navier Stokes equations (2.45)-(2.48). In the following, we look for the boundary 
layer {y < —h + 5 or y > h — S). The previous Thin Layer Navier Stokes equations 
simplify and we obtain the equations of acoustic boundary layer [Ch64]. 

• First we suppose that some reference state with null velocity is given. It is a 
priori the air at usual atmospheric pressure po and usual temperature ; with 
the state low of perfect gas, the reference density po is given. Second we search a 
field p(t, y), u{t, x, y), T{t^ x, y), p{t, x, y) of the form 

|p(t, X, y) = pq -\- p'{t, X, y) 

u{t, x,y) = + u'{t, X, y) 

T{t, X, y) = 6^ + T'{t, X, y) 

p{t, X, y) = pq + p'{t, X, y) . 

We linearize the equations (2.45)-(2.48) around the reference state (po, Po? ^o)- 
We suppose that we can neglect the nonlinear contributions of the bondary layer 
for studing the equations (3.21)-(3.24) of the main fiow. We recall that these equa- 
tions have been obtained by integrating the Thin Layer Navier Stokes equations 
on the complete width of the channel and the detailed analysis of the different con- 
tributions have been derived by Msallam [Ms98]. Recall that by doing this linear 
approximation, we neglect the acoustic streaming effect (see e.g. Batchelor [Ba67], 
Makarov and Ochmann [Mo97]), all separation effects inside the boundary layer 
(Merkli and Thoman [MT75]) and all random unstationary effects of turbulence 
[MT75]. The algebra is classical (Chester [Ch64]) and straightforward. We obtain 

v' p' T' 

Po Po Oo 

do' 

(4.3) + divu' = 0. 

du' d^u' dp' 

(4.4) po^-^— = - — 
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,^ r. dT' , d^r dp' 

^'■'^ = 

We note (with Chester [Ch64]) that equations (4.4) and (4.5) are two heat equa- 
tions coupled via the right hand sides. 

• The boundary conditions associated to (say) the bottom of the boundary layer 
{y = —h) are : 

(4.6) w'(t, x,y= -h) = 0. 

(4.7) r{t,x,y = -h) = 0. 

At the top of the boundary layer [y « —h -\- 5), we must mutch the boundary 
layer flow with the main flow : 

(4.8) w'(t, X, y — y —h + 5) — > (velocity in the main flow)(t, x) 

(4.9) T'{t, rr, y — > —h-\- 5) + — > (temperature in the main flow)(t, x) . 



5) The coupled problem. 

• We couple in this section the main flow in the pipe described in section 3 with 
the acoustic boundary layer presented in part 4. More precisely, the main flow is 
described by three unknown functions (density, velocity, internal energy) : 

(5.1) [0, +oo[x[0, L] 3 (t, x) I — y (p{t,x), u{t,x), e{t,x)) e [0, +oo[xIRx[0, +oo[ 

which represent the mean value in the section of the pipe of each field (denoted 
with a tilda in section 3). In the bounday layer, we suppose that the faces y = ±/i 
are composed by symmetric flows and we fix some transverse variable 7] e [0, +00 [. 





u(t,x) 








► 







Figure 2 Velocity field u{x, t) in the mean flow 
and velocity field ^(t, x, rj) inside the boundary layer. 

The unknowns are velocity ^ and temperature 6 in the boundary layer : 

(5.2) [0, +(X)[x [0, L] X [0, +oo[9 (t, x, rj) 1 — > (^(t, x, rj), e{t, x, r/)) e IR x [0, +cx)[ 
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Notice the important point concerning the modehing : we consider on one hand 
two velocity fields u and ^ and on the other hand two temperature fields T = 
e/Cy and 6. 

• The transverse scale 7] for describing the boundary layer fiow is very small 
compared to the transverse dimension h of the flow. Then it is consistent to set 
boundary conditions for 77 — y +00 : 



(5.3) 



— (t, X, rf) — )■ when ri — )■ +00 
orj 

— (t, X, rj) — > when 77 — > +00 



drj 

and we have done this particular choice in our simulations. Nevertheless, stronger 
boundary conditions for rj — )■ +00 that are compatible with the observed solu- 
tions in our numerical experiments could be the following ones : 

J <^(t, a;, if) — > u{t, x) when rj — y +00 
]^ 6{t, X, rj) — y T{t, x) when rj — > +00 . 



(5.4) 



For ?7 = 0, we just have to consider Dirichlet boundary conditions 

/ C(t, x,0) = 

\ eit, X, 0) = ^0 

where ^0 is the value of imposed temperature on the walls. 



(5.5) 



• The partial differential equations for the evolution of main flow variables are 
simply derived from equations (3.21)-(3.24) ; the source terms of stress and thermal 
flux at the wall in the right hand side of equations (3.23) and (3.24) are nomore 
obtained by solving the Thin Layer Navier Stokes equations (2.45)-(2.48) but the 
ones coming from the boundary layer model (4. 2)- (4. 5). We obtain in this way : 

(5.6) p(t,a:) = (7-l)pe 

d f ^ i\ d ^ 1-^ \ k 89 ^ . 

(5.9) —(pe + -pu ) + -Q^ypue + -pu + pu) = --— (t, x, 0) . 

• The evolution equations for boundary layer variables are obtained in a similar 
way from the heat equations (4.4)-(4.5) that model an acoustic boundary layer. 
With our coupled model, the pressure term comes from the main one- dimensional 
model and no more from the Thin Layer Navier Stokes equations ; it is therefore 
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considered as a source term and for this reason is placed on the right hand side of 
the equations. We get 

(5.10) 
(5.11) 







^ drf 


dp 
dx 


Po Cp 


89 

'dt 


drf' 


dp 

'dt 



• We observe that equation (5.10) is a dynamic equation that allow the predic- 
tion of velocity field (^(t, x, rf) as long as pressure field p(t, x) is known. It is not 
so clear for temperature equation (5.11) for the variable 9 due to the dymamic 
term ^ on the right hand side. Nevertheless, following a remark of Brenier 
[Br97], system (5.5) (second equation) and (5.11) can be replaced by the new 
unknow function a 

(5.12) a = poCpO - p 

that satisfies the following heat equation 



da k d 



(T 



dt Po Cp dr]"^ 

dp 

due to the fact that — = 0. The following nonhomogeneous Dirichlet boundary 

dr] 

condition is valid at the bottom of the boundary layer : 
(5.14) a{t, X, 0) = PqCpOq - p{t,x) . 

• It is clear that on one hand, that stress viscous term ^ || (r/ = 0) and thermal 
flux k ^{rj = 0) forces the main flow equations (5.7)-(5.9) and on the other hand 
that pressure fleld in the inviscid flow forces the boundary layer equations (5.10)- 
(5.11). We insist again on the fact that the main originality of our coupled model 
(5.6)-(5.11) consists of choosing two independent unknowns functions for velocity 
(in the main flow and in the bounday layer) and also two independent unknowns 
functions for temperature. We do not make the tentative to determine explicitely 
the boundary layer thickness 5{t^x) or the displacement thickness S*{tjX) (see 
e.g. Le Balleur [LB80]) in the way we set the coupled problem. In our approach, 
the boundary layer thickness for momentum and energy can be evaluated as a 
global (and nontrivial) result from the entire knowledge of functions [0, -|-oo[9 
r] I — y ^(t, X, f]) and [0, -|-oo[ 3 r] i — y 9{t, x, ry). 



• We recall briefly also the inflow-outflow boundary conditions at a: = and 
X = L concerning the mean flow variables. At the inflow {x = 0), the flow is 
subsonic then two conditions have to be considered : for axample, we give on one 
hand some data concerning the input velocity fleld uo{t) or the input pressure 
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field 7ro(t) and on the other hand the fact that entropy is not dissipated at the 
entrance of the channel : 

(5.15) u(t,0) = uo{t) or p{t , 0) = no{t) 

At the outflow, the mean fleld remains subsonic and the theory of characteristics 
(see e.g. Kreiss [Kr70]) show that only one scalar boundary condition is sufficient 
to set correctly the problem ; we choice nonreffecting boundary conditions (see 
Whitham [Wh74] or Hedstrom [He79]) : the outgoing wave is a so-called C+ 

P 

simple wave [Wh74] i.e. both specific entropy S = — and Riemann invariant 

2c ^ 

R = u take constant values everywhere in this part of the fiow. In 

7-1 

consequence this Riemann invariant R is advected with all characteristic celerities 
without distorsion and in particular the one with u — c velocity : 



(5.17) 




6) Generalization to axisymmetric geometry. 

• In this section, the pipe is nomore a two-dimensional channel but a three- 
dimensional cylinder with an axisymmetric geometry. The length of the pipe is 
still denoted by L and the letter h is used for the radius instead of half of the 

section. Hypothesis (2.1) concerning ratio — remains valid and we have 
(6.1) ^ « 1- 

As is section 2, Thin Layer Navier Stokes equations are a good approximation 
of the flow inside the entire geometry and this model takes now the following 
algebraic form in this axisymmetric geometry : 

at ox y oy^ ' 

. d . . 9(1 \ 1^/ \ p d f du\ 
(6-3) -^\pu) + -^[pu + p) + -^(puvy) = -^[y^] 
ot ox^ ' y oy y oy \ oy J 



6.4 - (^py-^y + py) = ^ — 

y oy^ y dy \ 



Idu 4y dv^ 
S'dx ^ 'S"dy^ 
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(6.5) 
(6.6) 



d f f 1 . 



+ 



d_ 

dx 



+ 

< y < h 



k d / dT 
y dy\ dy J y d y \ d y 



LL d f du 
u 



• The derivation of the coupled model can be conducted as in the previous 
sections. We first introduce the mean value of density, momentum and internal 
energy as in (3.2), (3.3) and (3.4) : 

(6.7) p{t,x) = — p{t,x,y)ydy 



(6.8) p{t,x)u{t,x) = ^ {pu){t,x,y)ydy 

(6.9) p{t,x)e{t,x) + -p{t,x)u^{t,x) = j (pe + ^pu'^)it,x,y) y dy 



We multiply equation (6.2) by y, integrate between and h and divide by /i^ : 2 
We get 

r -\y = h 
pvy = 0. 

L J?/ = 



dp d 2 



dt ' dx^^"^ ' 

The term inside the brakets in the left hand side of the previous relation is null 
due to the no slip boundary condition : 

(6.10) u{t, X, h) = v{t, X, h) = 0, 

and the conservation of mass becomes 



(6.11) 



dp d . — . 
+ —(pu) 



dt 



dx 







We make the same operation for the momentum equation (6.3) 



(6.12) < 



+ 



2 




y = h 


2ii 


" dw 


y 


= h 




pyuv 


y = o 




L dy\ 


y 


= 



Under the same hypotheses concerning the boundary layer presented in section 3, 
we have : 

2 

(6.13) (7-l)pe ~ — p{t,x,y)ydy 
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2 

(6.14) pu'^ c:i — j (pu^) {t,x,y)ydy 



We insert these evaluations inside relation (6.12) and obtain 
with 

(6.16) pit,x) = (7-l)pe. 

The treatment of the energy equation (6.5) is obtained by the same way, due to 
the boundary condition (6.10), approximations (6.13), (6.14) and 

X_ 2 / \ 

(6.17) pu(e + -u^) + pu c:i — ypu[e+ —) +puj{t,x,y)ydy 

(6.18) -(pe+-pu-) + -[pu(e + -u-)+VU^ = y^(M.ft). 

• We change the notations, replace the "tilde" unknown functions by letters 
without tilda, denote by ^(t, x, rj) the velocity in the boundary layer (0 < 77 < 
+00) and by 9{t, x, rf) the temperature in the same conditions. Due to the relation 

(6.19) y = h — T] , rj ^ 

we have to change the sign in the right hand side of equations (6.15) and (6.18). 
We get finally, as in (5.6)- (5.9) : 

(6.20) p{t,x) = (7-l)pe 

If + = ° 

d 1 d 1 2k 89 

(6.23) —{pe+-pu') + —(pue+-pu'+pu) = --—{t,x,0). 

We remark that there is just a factor of 2 that makes different the set of equatons 
(6.22) (6.23) from the set of relations (5.8) (5.9). 

• Inside the boundary layer, section 4 can be applied in a straightforward man- 
ner. We denote by u' the (infinitesimal) velocity, by p' the difference between 
pressure field p and ambiant pressure po and by T' the difference T — 9q. We 
have, due to the relations (4.4) and (4.5) : 

(6.24) p^'^-^-Uy'-f) = -'^ 

at y oy\ oy J ox 
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dT' k d ( dT'\ dp' 



(6.25) p,C,— ---[y—) 



y dy\ dy J dt 
In the boundary layer, we have y h and the curvature effects due to geometry 
are associated to the radius h and this distance is very big compared with the 
thickness 5 of the boundary layer : 

(6.26) 5 « h. 

To fix the ideas, velocity field u' can be expanded with an ansatz of the type 

(6.27) u' = Uf{y-) 

where /(•) is a regular function satisfying the conditions 

(6.28) /'(O) « r(0) « /(I) « 0(1). 

1 du' U 1 d^u' U 1 

Then t ~ -r^f'{^)i -^tt ~ -rT/"(0) and due to the relation 
h oy ho oy^ do 

1 du' 

(6.26), the term — — — can be neglected in comparison with the second term 

h oy 

d^u' 



dy'' 



• Finally the equations in the boundary layer can be written as 

(6.29) ^ - = -g^ 

(6.30) p,C,- -k^ = - 

as in the two dimensional case. The coupled problem in the axisymmetric case is 
composed by the set of equations (6.20)- (6.23) and (6.29)- (6.30). 



7) Numerical approximation of the coupled problem. 

• The coupled system defined in section 5 is composed by five partial differential 
equations (5.7)-(5.11), the state low of perfect gas (5.6), the boundary conditions 
(5.3) (5.4) at the top-bottom of the pipe and by the inflow-outflow boundary condi- 
tions (5.15)-(5.17). We discretize this system of equations in the following manner. 

• First we introduce some integer J and an associated space step Ax : 
(7.1) = J 

and some time step At is chosen below. We define the discrete variables pj, uj, ej 
for density, velocity and energy at discrete point xj = j Ax {j = 0, 1, 2, . . . , J) 
and at time = mAt : 
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(7.5) G{W) = -(^0,f3f,-^(t,x,0),(3k — (t,x,0)) 



pj" p (m At , j Ax) 

(7.2) { uj ^ u{mAt, j Ax) 
e (m At , j Ax) 

and we suppose that state equation is satisfied at time step mAt and at vertex 
Xj = j Ax : 

(7.3) = (7-l)p7e7, 0<j<J, 0<m<n. 

• We introduce the conservative variables W for mean flow : 

p , pu , p[e + —)j , 
the physical flux function f{W) : 

pu, pu^ +p , pu(e+—) +puj 
and the source term due to the boundary layer : 

^(t, X, 0) , (3k^ 
or] ' or] 

where the variable (3 is defined by the condition 

(7 fi\ a — j ^ the plane case of relations (5.7)-(5.9) 

^ ' ^ I 2 in the axisymmetric case of relations (6.21)-(6.23). 

Then the equations (5.7)-(5.9) and (6.21)-(6.23) can be written in a more compact 
form : 

dW d 

• Variables (7.2) (for j = 1, 2, . . . , J — 1) are advanced between times = 
nAt and t""^^ = (n + l)At according to the Lax-Wendrofi^ [LW60] numerical 
scheme. This scheme is founded on a second order Taylor expansion in time of the 
conserved variables : 

(.3) -r--r.-(?)>^-^(5); 

that is exact up to a third order term relatively to variable At which is omitted 
in the numerical scheme (7.8). The first derivative in time I 1 is directly 



dt 



evaluated thanks to equation (7.7) 

'9WY _ (q^^^v (am)\ 



3 



and more precisely 



Regis Msallam and FRANgois Dubois 



n 



= < 



The second derivative in time 



/ df df) \t 

(O , ^ (r, j Ax, 0),/3k — (r, j Ax, 0)) 

(^pu , pv? + p, pu{e+ ^) + puj 
[pu , py?' + p, pu{e+ ^) + pu^ 



2 Ax 



is obtained by a derivation of equation 



(7.7) that takes into account the Schwarz property for partial derivatives 



dtG{w)y - 



l_ ( df{w) \ 

dx\ dt ) 



d 



and after discretization of the — operator by finite differences, we get 

ox 



(7.10) 



dtG{W) 



n 



- K/(w')),_,/,-(a*w);-v2 



We use classical expressions for the discrete operators presented in equation (7.10) : 
the time derivative of right hand side of equation (7.7) is local in space and will 
be evaluated "more above" : 



(7.11) [dtG{W))\ 



I (O , /3/. I (r, , A.:, 0) , /^fc g (r, , A.:, O)) 



t 



the jacobian matrix ypw fiy^)j ^ at the intermediate point (j + l/2)Aa: is 
evaluated thanks to a simple two-point mean value formula : 



(7.12) {dwf{W) 



J+l/2 



1 

2 

+ 



d[pu, pv? + p, pu{e+ ^) +pu)^^ " 
d(p,pu,p{e+^))' 

2 

d[pu , pu^ + p, pu(e+ ^j-) + pu) 
d{p,pu,p{e+^)f 



+ 



3 



and the time derivative of conservative variables at intermediate point (j + l/2)Aa: 
is obtained with a centered scheme : 
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(7.13) (9tH^);+i/, 



dr] 



(O, /3/.^(r,0- + l)A:., 0),/3/c-(r, 0- + l)A:., O)) 
ypu , pu'^ + p, pu{e+—) + puj j 



Ax 



i+i 



• The source term G{W) (relation (7.5)) is simple to represent with an integral 
formula, due to the simple structure of heat equations (5.10) and (5.11). We have 
(see e.g. Morse and Feshbach [MF53]) : 

(7.14) erf(x) = ^ / e"^ da , 



(7.15) ^{t,x,rj) = -- ^{z,x)erf 

Po JQ ox 



po 



(7.16) 0{t,x,r]) = $0 + 



dz 



7) 



k 



PoCj 



it-z) 



dz , 



and after derivation relatively to transverse variable r] 



(7.17) §j(t,x,0) 



p Jo dx^^ ^' \l 



P 



dz 



4 



dp 



pOTTZ 



P 



dz 



do 1 
(7.18) — (t,x,0) = -/ ^{t-z,x) 

or] k jQ dt y poCpTTZ 

In consequence, the source term (^G{W)^ is numerically evaluated according to 



/ 



(7.19) {GiW)y_ = 







n — 1 



m=0 



\ 



/.(m+1) At 


dp 


JmAt 


dx 


/.(m+1) At 


dp 
dV 


JmAt 



\ 



(n At — z, Xj) ^ 



P 



dz 



pOTTZ 



P 



„ dz , 

PoCpTTZ ) 
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The intermediate integrals in the second hne of right hand side of (7.19) is ap- 

dz 

proached with a two-point quadrature formula relatively to the measure 



(7.20) 



'a 



dz 1 



dz 



'a 



(ip{a) + ip{b)) 



b — a 



and integrals in third line of relation (7.19) are numerically approached by a one- 
point quadrature formula : 



(7.21) 



(p{z) 



'a 



dz 

7^ 



fa + b 



'a 



dz ^ f a + b\ b — a 



V 2 / ■ 

We deduce, due to quadrature relations (7.20)-(7.21) and elementary use of finite 
differences : 



f (G2(W'))J = 



m=0 



2Ax 



G:i{W) 



and finally 



/ n--m-l I n--m\ _ / n--m-l , n--m\ 
1 



J- 
n-1 



1 



2/3 E 



m=0 



P, 



P 



n—m—1 
3 



Po ^/m -\- \Jm-\-\ 
1 



(7.22) {g{W)Y = fo, {G2{W)y , [Gs{W)J 




x=0 



Figure 3 Characteristic directions at the entrance x = 0. 

From previous evaluations, the time derivative of the source term (^dtG{W))^ is 
computed with a simple first order scheme : 
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(7.23) mW))" = (O , - "^^'^"'^ - 



At ' At y 

• We neglect the boundary layers when considering numerically the boundary 
conditions at the input and at the output of the domain. The boundary conditions 
at j = and j = J are numerically implemented using the method of charac- 
teristics (see e.g. Whitham [Wh74]). We distinguish between threen cases : input 
pressure wave, input simple velocity wave and nonrefiecting output. In the case of 
an input pressure wave (at j = 0), two characteristics directions are going inside 
the computational domain (for celerities u and w + c) and one (associated with 
celerity w — c) is going outside (see Figure 3). We wish to define the state W,^"'"''" 
at the first mesh point and at time n + 1 ; all the states at time level n are 
supposed to be given and the pressure field at time level n + 1 is imposed to be 
equal to some numerical value tt""*"^ due to the boundary condition. We denote 
by Co, Po and respectively the sound celerity, the pressure and the entropy of 
the air at rest at usual conditions of temperature and pressure. We first determine 
an external sound celerity Cg and an external velocity Ue associated with a C+ 
input wave ; we have classically from locally linearized theory [Wh74] : 

(7.24) «e = ^^^^^ 

po Co 

7-1 

(7.25) Ce = Co H — Ue . 

Secondly we interpolate data at time level n and at the foot-point P going back- 
ward along the u — c characteristics starting at t'^'^^ from x = : 

(7.26) Wp = (l - ^^°^) + . 

The state VF(5I^''"^ is finally defined by the following three conditions : the charac- 
teristic variable associated to the u — c wave is constant between states Wp and 

(7 07\ -,,"+1 ^^0^^ — „ "^^P 

Uq - Up -, 

7 — 1 7 — 1 

the entropy of state Wp is equal to the entropy at rest : 

(7.28) S^+' = So 

and the characteristic variable associated to the u -\- c wave is constant between 
external state and state W,^"*"^ : 

(7.29) + ^ = „^ + . 

7 — 1 7 — 1 
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With this implementation, the single pressure datum variable 7r"'+^ is used for 
two incoming waves and the outgoing wave is not reflected. We remark that, as in 
[DF89], nothing in what we have done imposes strongly the condition p{Wq^^) = 
7r'^+^. In some sense, this boundary condition is transparent to the outgoing waves. 



• We use the same notations for the input simple velocity wave associated to 
datum at the time level under study. This datum is supposed to be suffi- 

ciently small in order to be considered as an acoustic velocity. We first determine 
the celerity of an external state by a relation similar to (7.24) : 

(7.30) Ce = co + ^f/"+S 

we interpolate a state Wp at the foot of the u — c characteristic direction using 
relation (7.26) and the boundary state Wq^^ is computed according to relations 
(7.27) along the outgoing characteristic, (7.28) along the u characteristic direction 
and the following relation along the w + c incoming characteristic : 

(7.31) + = ^""^ 



7 — 1 7 — 1 




Figure 4 Nonrefiecting output at x = L. 



• For a nonrefiecting output at x = L and j = J (see Figure 4), the external state 
is the air at rest and is obtained by going backward along the u — c characteristic 
direction : 



(7.32) u 



_ 2c;+^ _ 2co 
~ 7-1 ~ ~7 - 1 



n+l ^"-J 



Concerning the waves going outside the computational domain, we define the foot- 
point Q associated with the u-\-c characteristic direction with the same idea than 
previously : 
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(7.33) Wq = (i- wu + wy 

and we say that the associated characteristic variable is constant between this 
state Wq and state Wj"^^ : 

(7.34) u-/' + = + 



7 — 1 7 — 1 

We suppose also than a relation similar to (7.28) determines the entropy at the 
limiting vertex : 

(7.35) 5^+^ = 50. 



8) First test cases. 

8.1) Nonlinear perfect oscillating fluid. 

• In this sub-section, we neglect all the viscous effects. The continuous model 
is given by equations (5.6)- (5.9) with fi = and k = and all the discrete 
equations correspond to Lax-Wendroff scheme (7.8)- (7. 13) without source terms. 
The first test case consists of a simple wave going inside the domain {a; > 0}. At 
time t equal to zero, the fluid is at rest (with pressure po and temperature To 
that correspond to usual thermodynamics conditions) and at x = 0, a source of 
velocity w(0, t) is supposed to be given. It defines a simple wave, submitted to 
hypothesis 

8.1) u = - 

7 — 1 7 — 1 

and if a; = X{t) is the solution of the differential equation that defines the char- 
acteristic line, i.e. 

(8.2) ^ = « + c, 

we have (see e.g. Whitham [Wh74]) 
2 c 

(8.3) u + = Cste. 

7-1 

By elimination of sound celerity c between equations (8.1) and (8.3), velocity u{») 
has a constant value along the characteristic (8.2)- (8.3), and it is also the case for 
sound celerity. We deduce that u-\- c depends only of its value for x = and the 
slope of characteristic direction is constant : 

(8.4) u + c = Co + ^^wo(to). 

Then characteristic lines are straight lines. Moreover, if a: = L is some given 
abscissa, the time II — to for the wave to propagate between a: = at time 
t = to and a: = L at time t = tL is given according to the following relation : 
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(8.5) 



tL — to 



I. 



u + c 



Co H ^^o(^o) 



• We construct velocity field at the particular station x = L with regular time 
steps r. We search velocity w(L, nr) according to the relation : 

(8.6) m(L, nr) = wo(to,n) 

where to, n is the solution of the following nonlinear equation : 

L 



(8.7) nr - to, 



Co H 7;—uo{to,n) 



Equation (8.7) is solved with the help of a Newton algorithm detailed in [Ms98] 
for a sinusoidal input velocity Uo{0) : 

(8.8) ^0(6*) = Uosm{uJoe) 

and the Newton algorithm is congergent without any problem as long as the char- 
acteristic lines does not intersect, i.e. under the condition 

2 Co 

(8.9) L < I/ghock ^ 



(7 + 1) ujQ Uo 

We introduce adimensionalized abscissa s relatively to I/ghock 
(8.10) s ^ 



-Ishock 



The output velocity u{L, •) is a periodic function of time with period T = 27r/uJo 
and parameter r has been chosen such that 

with a big integer N (of the order of 10 typically) in order to proceed a precise 
signal treatment. The nonlinear distorsion effect induces harmonics kujo {k = 
2, 3, • • •) of fondamental pulsation uq and they are predicted up to = 40. Note 
that the time step r for computing exact solution has been chosen sufficiently 
small in order to avoid aliasing effects when computing the fast Fourier transform. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
adimensionalized time t/To 



Figure 5 Signal at abscissa s = 0.8 
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Figure 6 Signal for three abscissae 
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Figure 7 Magnitude spectrum for three abscissae 




Figure 8 Influence of the space step Aa: 



• We first test the effect of non-absorbing boundary condition on the numerical 
flow computed with help of La:x-Wendroff scheme inside the domain. We consider 
two simulations on two computational domains with the same space step Aa:. 
One domain is of lenght L and the other one is constructed in order to be sure 
that the boundary scheme is not active at x = L. Then we check that nonlinear 
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treatment (7.32)- (7.35) induces, for waves that compose the distorted signal at 
X = relative errors in velocity that are inferior to 2 % for waves containing 
more than 10 grid points. 

• We compare this simple wave with the numerical solution computed with 
Lax-Wendroff scheme. We choose a simple sinusoidal velocity profile (8.8) at the 
inflow (see also (7.30) and (7.31) for the complementary boundary conditions at 
the inflow) and a non-reflecting boundary condition at x = L (see also relations 
(7.31) to (7.35)). On Figure 5, we plot temporal output signal for exact (charac- 
teristics) and approached (Lax Wendroff) methods at station s = 0.8. We notice 
that Lax-Wendroff scheme is correct for prediction of this kind of nonlinear wave. 
We recover the distorsion of the wave with a proflle more and more sharp as vari- 
able s is increasing. We compare three results for the Lax Wendroff scheme at 
s = 0.005, 0.4 and 0.8 on Figure 6 and the associated spectra for these three loca- 
tions on Figure 7. We verify on Figure 7 that distorsion induces an enrichment of 
spectrum with a transfer of energy from low frequency to higher frequencies. More- 
over, comparison of spectra for both methods shows that Lax-Wendroff scheme is 
operational for good prediction of output signal. We compare also spectra of out- 
put signals for different values of space steps Ax with constant CFL number that 
induces proportional values of At. For this particular simulation (wo(«) given by 
relation (8.8) and s = 0.8) we observe (Figure 8) that numerical damping is com- 
patible with harmonic k = lb for 100 temporal points by period {i.e. j| ?a 6.3 
points for this particular harmonic) and with harmonic = 25 for 190 points by 
time period (^ = 7.6 points for one period). 

8.2) Linear wave with visco-t hernial boundary layer effects. 

• In this sub-section, we compare our numerical model with the linear Kirchhoff 
theory obtained by linearizing convective effects around a null velocity. We refer 
to Bruneau et al [BHKP89] for this classical approach in the context of first order 
theory with thin boundary layer. Recall that a wave with pulsation u can be a 
particular solution of linear Kirchhoff theory if the phase 

(8.12) $ = ut - Kx 

admits a dispersion relation of the type 



(8.13) 



K = 




i a{uj) 



with 



(8.14) 
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(8.15) aico) = (.[^ + (7 - 1) J^^) ^ . 

As long as the wave propagate, there is dispersion and damping of this wave. 
Dispersion is due to the fact that local phase velocity c'{u3) depends on frequency 
(see relation (8.14)). Damping is associated to the real part of the constant of 
propagation and a{u)) is a damping coefficient. 



45 




number of points per wavelength (log scale) 



Figure 9 Relative error on amplitude between 
the Kirchhoff theory and the Lax-Wendroff code 



• Sinusoidal wave of small amplitude have been simulated by the numerical 
model. In order to avoid the essential of nonlinear effects, a very small amplitude 
has been chosen for the wave. We simulate the propagation of an acoustic wave. 
The sinusoidal input profile is computed over a distance L physically of the order 
of 1 meter in a pipe of diameter of the order of 1 centimeter. We compare both 
amplitudes of the waves in the first case by using Kirchhoff model and in the 
second case with the pure numerical solver. The relative errors of the predicted 
amplitude are plotted on Figure 9. We use as space variable the number of grid 
points for one lenght wave. With more than 25 points by lenght wave, we observe 
that relative error for pressure field is inferior to 5%. The results concerning the 
phase obay to the same conclusion. 



Model for coupling a perfect flow with an acoustic boundary layer 




Figure 10 Comparison of flow field at s = 0.8 with and without loses 



8.3) Combined nonlinear propagation and linear boundary layer. 

• In this sub-section, we compare the shape of the wave with and without visco- 
thermal boundary layer. Observe that view results are available in the literature 
for this kind of elementary coupled problem. First comparison have been done 
with results obtained independently by Menguy and Gilbert [MG97b]. Figure 
10 presents a temporal signal of velocity at fixed abscissa s = 0.8. There is an 
important damping of the wave and we recover that wavefront is less sharp with the 
presence of the boundary layer. Notice here the important remark that viscosity 
associated to the boundary layer is much more important that the one due to the 
thin layer Navier Stokes equations. 

8.4) Trombone modelling. 

• In [HGMW96], Hirschberg, Gilbert, Wijnands and the first author have demon- 
strated experimentally that for high level of amplitude {forte, fortissimo) , there 
are important nonlinear propagation effects in the trombone which can lead to 
shock waves. From a modelling point of view, the slide of a trombone can be 
viewed in first approximation as unidimensional pipe of lenght of the order of 1.5 
meter and redius h = 7 mm. A typical incident pressure wave at the entrance 
of the slide is propagating along the slide. It is a low frequency signal. At the 
output, we use an absorbing boundary condition. In fact, a complete model of 
trombone would include the discretization of the bell. But in the strong fiairing 
part of the bell, the flow is no more quasi-unidimensional and our model is no more 
relevant (see Amir, Pagneux and Kergomard [APK97]). Because the slide is the 
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largest lenght with a cylindrical shape of the instrument, we conjecture that the 
essential of nonlinear effects occur in this part of the instrument. We restrict our 
simulations to the slide alone and make the hypothesis that nonlinear interaction 
between incident and reflected waves are negligeable. This hypothesis has been 
verified with numerical tests [Ms98] . 
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Figure 11 Magnitude spectrum at the output of the slide 

• We have synthetiscd a typical input signal containing four harmonics. This 
signal has been propagated with a linear propagation with losses and with non- 
linear advection with and without losses. Figure 11 represents these three output 
signals. There is of course no creation of superior harmonics with linear propaga- 
tion (done with linearized computer software developed by Quinnez [Qu95]). We 
recover the four initial input modes and damping is visible (2 dB) on the sound 
pressure level in Figure 11. With nonlinear propagation without losses, new har- 
monics for k > b are created. Moreover the amplitude of all the harmonics (except 
the first one) is amplified by nonlinear propagation. This corresponds to trans- 
fer of energy from lowest frequency to higher frequencies in order to go towards 
thermostatics equilibrium where we have equi-partition of the energy between all 
the modes. With both nonlinear and linear effects, previous results are damped 
with an amplitude varying between 1 and 5 dB for the 8 first modes. The effect 
of losses compensates the one of nonlinear propagation. Nevertheless, nonlinear 
effects dominate the dynamics. For example the amplitude of the 4th mode is 
increased of 5 dB compared to the input signal. 
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• Our numerical results confirm previous experiments done in Eindhoven : trom- 
bone's radiated sound is enriched by nonlinear effects which occurs in the instru- 
ment. This effect is known by the musicians as the "brassy" sound, typical for the 
trombone at loud tones (see e.g. [GM98]). Some sonor amplitude examples are 
available on the net at the following http://www.icp.inpg.fr/~pelorson/sons.html. 



9) Conclusion, acknowlegments. 

• In this study, we have proposed to use the Thin Layer Navier Stokes equations 
as primitive ones to study propagation effects in thin pipes. This complete model 
neglects diffusive effects in the stream direction. Second, we have derived from this 
primitive set of partial differential equations a coupled model of five equations that 
takes into account both nonlinear uni-dimensional propagation and linear diffusion 
in acoustic linear boundary layer. The numerical coupling of this two models have 
been done and the approach is original : there is no explicit need of the displace- 
ment thickness but a set of two velocities and two temperatures (one in the main 
nonviscous flow and one in the boundary layer) allow this coupling. First numer- 
ical experiments have shown global coherence with previous classical models in 
computational acoustics (characteristics, linear Kirchhoff theory). Moreover, first 
application to one-dimensional modelling of trombone confirm the importance of 
nonlinear wave propagation and in particular the "brassy" sound that is famil- 
iar to jazz musicians. The extensions of this work concern a complete treatment 
of nonlinear waves with precise simulation of shock waves in the trombone, new 
numerical experiments where convolution effects in the boundary layer are com- 
puted with direct numerical resolution of heat equation, coupling for aerodynamic 
flows where displacement effects play an important role (see e.g. Lagree [La2k]), 
and mathematical study of simplified models. The authors thank P.Y. Lagree for 
precise reading and comments on the first version of the manuscript. 
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